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Abstract 

The method recently developed to include Van der Waals interactions in the Density Functional 
Theory by using the Maximally-Localized Wannier functions, is improved and extended to the case 
of atoms and fragments weakly bonded (physisorbed) to metal and semimetal surfaces, thus opening 
the way to realistic simulations of surface-physics processes, where Van der Waals interactions play 
a key role. Successful applications to the case of Ar on graphite and on the Al(100) surface, and 
of the H2 molecule on Al(100) are presented. 
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Understanding adsorption processes is essential to design and optimize countless material 
applications, and to interpret scattering experiments and atomic-force microscopy. In par- 
ticular, the adsorption of rare-gas atoms on metal and semi-metal surfaces is prototypical [l| 
for physisorption. The weak binding of physisorbed closed electron-shell atoms or saturated 
molecules like H 2 is due to an equilibrium between attractive long-range Van der Waals 
(VdW) interactions and short-range Pauli repulsion. Density Functional Theory (DFT) is 
a well-established computational approach to study the structural and electronic properties 
of condensed matter systems from first principles, and, in particular, to elucidate complex 
surface processes such as adsorptions, catalytic reactions, and diffusive motions. Although 
current density functionals are able to describe quantitatively several systems at much lower 
computational cost than other first principles methods, they fail to do so[2j for the descrip- 
tion of VdW interactions, particularly the leading R~ 6 term originating from correlated 
instantaneous dipole fluctuations. The key issue is finding an accurate way to include VdW 
effects in DFT without dramatically increasing the computational cost. We have recently 
proposed 3] a novel approach, based on the use of Wannier functions, that combines the sim- 
plicity of a semiempirical formalism [4j with the accuracy of the first principles approaches 
(see for instance ref.|5(), and appears to be promising, being simple, efficient, accurate, and 
transferable (charge polarization effects are naturally included). The results of test appli- 



cations to small molecules and bulk graphite were successful [3]. In this paper we describe 
improvements in our method which allow to extend it to the case of the interaction between 
an atom or a neutral fragment and a surface. 

Crucial to our analysis is the use of the Maximally-Localized Wannier function (MLWF) 
formalism [fl, that allows the total electronic density to be partitioned, in a chemically trans- 
parent and unambiguous way, into individual fragment contributions, even in periodically- 
repeated systems. The MLWFs, {w n (r)}, are generated by performing a unitary transfor- 
mation in the subspace of the occupied Kohn-Sham orbitals, obtained by a standard DFT 
calculation, so as to minimize the functional Q, defined as : 

n = Y. S n = (( w n\r 2 \w n ) - (w n \r\w n ) 2 ^j . (1) 

n n 

Besides its spread, S n , each MLWF is characterized also by its Wannier-function center 
(WFC); for instance, if periodic boundary conditions are used with a cubic supercell of side 
L, the coordinate x„ of the n-th WFC is defined roll as 



x n = -7^ Im ln(wn|e %2 Z x \Wn) • (2) 

If spin degeneracy is exploited, every MLWF corresponds to 2 paired electrons. Suitable 
codes[z] are available, which allow the efficient generation of the MLWFs, by adopting a 
proper Appoint sampling of the Brillouin Zone (BZ), which is crucial for metallic systems. 

Starting from these MLWFs the leading R~ 6 VdW correction term can be evaluated^ by 
making the reasonable assumption 8j of exponential localization of the MLWFs in real space, 
so that each of them is supposed to be an H-like, normalized function, centered around its 
WFC position, r n , with a spread S n : 

^(|r-r„|) = -^e-^ |r - rn| . (3) 

Then the binding energy of a system composed of two fragments is given by £& = Eo+Ey^w, 
where Eq is the binding energy obtained from a standard DFT calculation, while the VdW 
correction is assumed to have the form: 

^VdW = ~J2fnl ( r nl) , (4) 

n,l r nl 

where r n \ is the distance of the n-th WFC, of the first fragment, from the /-th WFC of the 
second one, the sum is over all the MLFWs of the two fragments, and the C§ n \ coefficients 
can be calculated directly from the basic information (center positions and spreads) given 
by the MLFWs. In fact, using the expression proposed by Andersson et al. (see Eq. (10) of 
ref-jsl) that describes the long-range interaction between two separated fragments of matter: 



C Gnl = _L_ [ dr ( dr' V /p " (r)p/(r/) = ^_ / dr [ dr > w n(r)Mr>) 

&n 327T 3 / 2 7|r|<r c J\r<\<r> c yf p n (r) + y/ Pl(r') 32tT 3 / 2 i|r|<r c J\r'\<r' c W n (r) + W 1 (t J ) 

(5) 

where p„(r) = u^(r) is the electronic density corresponding to the n-th MLWF, C 6n ; is given 
in a.u., and the r c , r' c cutoffs have been introduced P, [lO] to properly take into account both 
the limit of separated fragments and of distant disturbances in an electron gas. By using 
the analytic form (see Eq. (J3j)) of the MLWFs, it is straightforward [3] to obtain the cutoff 
expressed in terms of the MLWF spread: 
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S n y/3 (0.769 + 1/21x1(5;)) 



(6) 



and to evaluate very efficiently the multidimensional integral of Eq. ([5]); when each MLWF 
describes 2 paired electrons the Ce n i coefficients must be multiplied by a 72 factor 

In Eq. (BJ fniij) is a damping function to cutoff the unphysical behavior of the asymptotic 
VdW correction at small fragment separations. For this we have chosen the form[4j, |ll| : 



fm{r) 



(?) 



1 + exp(—a(r / R s — 1)) 

where[4] a ~ 20 (the results are almost independent on the particular value of this parame- 
ter), and R s = -Rvdw + -Rydw * s the sum °f the VdW radii of the MLWFs. In ref. |3( -Rvdw 
was determined as the radius of the 0.01 Bohr -3 electron density contour. For the present 
applications to extended systems with metal or semimetal surfaces, after extensive testing, 
we found that a better choice is to equate -Rvdw to the cutoff radius introduced in Eqs. (5) 
and (6), -Rydw — r c which has the additional advantage of not being dependent on an 
given electron density threshold. It should be stressed that the results reported in ref. 
relative to isolated fragments and bulk graphite, are essentially unchanged if recomputed by 
adopting this new -Rvdw definition. 

The Eq binding energy can be obtained from a standard DFT calculation (for instance, 



my 



using the Quantum-ESPRESSO 



proximation in the revPBE flavor 13| . This choice 



2] ab initio package), with the Generalized Gradient Ap- 



is motivated by the fact that revPBE 



is fitted to the exact Hartree-Fock exchange, so that the VdW correlation energy only comes 
from the VdW correction term, as described above, without any double-counting effect. The 
evaluation of the VdW correction as a post-DFT perturbation, using the revPBE electronic 
density distribution, represents an approximation because, in principle, a full self-consistent 
calculations should be performed; however recent investigations on different systems have 
shown that the effects due to the lack of self-consistency in VdW-corrected DFT schemes are 

n 

negligible. The method described above can be refined by considering the anisotropypl] of 
the MLWFs (details will be published elsewhere 15]); however, since previous calculations 
have shown that this has small effects, it has been not included in the present applications. 
Remarkably, the whole procedure of generating the MLWFs and evaluating the VdW correc- 
tions represents a negligible additional computational cost, compared to that of a standard 
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DFT calculation. 

We have applied our method to the case of Ar on graphite and on the Al(100) surface, 
and of the saturated H 2 molecule on Al(100). Adsorption on graphite has been modeled 
using an hexagonal supercell containing 36 C atom distributed over 2 graphene sheets, with 
a sampling of the BZ limited to the T point (preliminary tests and previous studies [la. \v\ 
show that these choices are adequate); in the case of the Al(100) surface the supercell was 
orthorhombic with a surface slab made of 32 Al atoms distributed over 4 layers, and a 2 x 2 x 1 
sampling of the BZ was used; no appreciable difference in the equilibrium properties was 
observed in test calculations with a thicker slab of 64 Al atoms over 8 layers (of course a 
thicker slab would be instead necessary to describe well the far-from-the-surface, asymptotic 
behavior, where the binding energy is expected to decay as z~ 3 , z being the fragment- surface 
distance). For a better accuracy, in these applications it has been necessary to modify our 
algorithm in such a way to include interactions of the MLWFs of the physisorbed fragments 
(Ar or H2) not only with the MLWFs of the underlying surface, within the reference supercell, 
but also with a sufficient number of periodically-repeated surface MLWFs (in any case, given 
the R~ G decay of the VdW interactions the convergence is rapidly achieved). Electron-ion 
interactions were described using norm-conserving pseudopotentials (in the case of Al only 
the 3 valence electrons per atom were explicitly included). In principle, for evaluating 
adsorption properties in periodically-repeated, asymmetric configurations, one should add 
a dipole correction 18j that compensates for the artificial dipole field introduced by the 
periodic boundary conditions; however we have checked that, in our cases, this correction is 
very small (just a few meV in the binding energy of Ar on Al(100)). 



The absorption o: 
sively over the years 



noble gases on graphite and on metal surfaces has been studied exten- 
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19] because it serves as the paradigm of weak adsorption. Actually, 



despite the conceived "simplicity" of these systems, even the most basic question (what is 
the preferred adsorption site ?) has not been answered in an entirely satisfactorily way. In 
principle, due to the non-directional character of the VdW interactions, sites that maximize 
the coordination of the adsorbate atom were expected, so that it was typically assumed 
that the adsorbate occupies the maximally coordinated hollow site. The actual scenario is 
more complex: for Xe and Kr a clear preference is found H| for adsorption on metallic 
surfaces in the low-coordination top sites (this behavior was attributed to the derealization 
of charge density that increases the repulsive effect at the hollow sites relative to the top site 
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and lifts the potential well upwards both in energy and height); for Ar the situation seems 
to be different: comparison 

the hollow sites are favored for Ar on Ag(lll) and on graphite, although, in this latter case, 
this configuration is preferred over two other possible sites (top and bridge), by only a few 



meV 



16|. 



Ar on Al(100) represents a critical test for our method, in fact the Al case is particularly 
challenging for a Wannier-based scheme since Al is the metal which most closely approxi- 



mates a free electron gas system: hence the electronic charge is relatively deloca^ 
assumption of exponential localization of the MLWFs is no longer strictly valid 



ized and the 
. However 



the following results show that, even in this case, our method works and this does not come 



to a surprise, 
also to metals 



n fact, on the one hand, the MLWF technique has been efficiently generalized 



21L 



221 ] , on the other, bonding in metallic clusters and in fee bulk metals (like 



2l|, 



Al) can be described in terms of H-like orbitals localized on tetrahedral interstitial sites 
which is just in line with the spirit of the present scheme. 

In the Tables I and II we report our computed binding energies and equilibrium fragment- 
surface distances, compared to the most reliable (to our knowledge) experimental and the- 
oretical reference data, and to the results of LDA calculations (we have also reported the 
values obtained in ref. 3| for Ar interacting with the benzene molecule). As can be seen, the 
general performance of the method is quite satisfactory; in fact, the improvement achieved 
by including the VdW correction, with respect to the pure revPBE scheme (which gives 
completely unphysical results, namely a potential well very small and located too far from 
the surface) is dramatic. 

In the case of Ar on graphite, the hollow configuration is energetically favored, although 
by just a few of meVs with respect to the other two configurations, in good agreement with 
previous studies[16|. Interestingly, our estimated Ar-graphite surface distance essentially 
coincides with the sum of the Ar and C VdW literature radii (1.88+1.73=3.61 A), a behavior 
experimentally observed in the related case of Xe adsorbed on graphite [19(. Moreover, the 
fact that the Ar-graphite distance is not appreciably smaller in the hollow site, compared to 
the top one, could be rationalized in terms of the potential lifting, due to increased repulsion, 
mentioned above. Note that the binding energy of Ar on graphite is considerably larger than 
that of the Ar-benzene complex, although the equilibrium distance is similar; this behavior 
is clearly due to the VdW interaction of Ar with the electronic charge outside the underlying 
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C ring, and is not reproduced by the LDA approach which favors short-range interactions. 

Concerning Ar on Al(100) (see also Fig. 1), specific experimental values are not available, 
however the experimental binding energy of Ar on several other metals is found to be in the 



range between 30 and 100 meV 



24| (between 70 and 85 meV|25j for noble metals) 



in agreement with our VdW-corrected results. We also mention old theoretical estimates of 
a binding energy of about 200 meV 26], and of 70 meV using a jellium model 271 ]. 



In t 



regime 



re case of H2 on Al(100) the molecule is essentially a free rotor in the physisorption 



281 ] and its interaction with the substrate exhibits only a slight anisotropy; moreover 



the effect of changing the position of the molecule with respect to the substrate is small, so 
that we report only the results relative to a single, representative configuration. Even for this 
extremely weakly bonded system the results are in good agreement with the reference values 
(we also mention that the binding energy of H2 on Mg is is predicted to be 17 mevj^]). 

Looking at the tables, on can see that the binding energies are reasonable reproduced by 
the LDA scheme, although this is actually accidental (the well-known LDA overbinding, due 
to the overestimate of the long-range part of the exchange contribution, somehow mimics the 
missing VdW interactions), moreover the equilibrium distances are clearly underestimated. 

In conclusion, we have extended our recently developed scheme, to include VdW interac- 
tions in the DFT by using the MLWFs, to the case of fragments weakly bonded (physisorbed) 
to metal and semimetal surfaces, and we have reported results of applications to the case 
of Ar on graphite and on the Al(100) surface, and of the H2 molecule on Al(100). The 
good performances of the method clearly indicate that it can be very useful to investigate 
many realistic surface-physics processes, where VdW interactions play a key role. Of a par- 



ating surfaces could be 
. Finally it must be 



ticular value is the possibility of dealing with metal surfaces (insu 
somehow treated even using atom-based semiempirical approaches 
stressed that a large area for future improvements of the method exists. In fact, different, 
more sophisticated schemes to utilize the MLWFs could be developed: for instance, one 
could adopt gaussians instead of exponential, H-like, functions, because multidimensional 
integrals are more easily evaluated; orbitals of symmetry different from the s-like one could 
be used for specific applications; partially occupied MLWFs 29( , with improved localization 
and symmetry properties, could be introduced (particularly for metallic systems); different 
damping functions, and improved, reference DFT functionals, with respect to revPBE, could 
be chosen,... 
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TABLE I: Binding energy, in meV, of Ar on graphite and Al(100) and H2 on Al(100), computed 
using the standard DFT-revPBE calculation, and including the VdW correction, compared to the 
LDA result, and available theoretical and experimental (in parenthesis) reference data. 



system 


revPBE 


revPBE+VdW 


LDA 


ref. 


Ar-gr. hollow 


-2 


-133 


-76 


-116 a , -111 6 (-119 c ) 


Ar-gr. top 


-2 


-131 


-69 


-106 6 


Ar-gr. bridge 


-2 


-128 


-69 


-107 6 


Ar-benzene 


-2 


-66 


-72 


-65 d (-49 e ) 


Ar-Al(lOO) hollow 


-3 


-72 


-71 




Ar-Al(lOO) top 


-3 


-71 


-66 




H 2 -A1(100) 


-2 


-20 


-24 


-19' (-28») 



"Reference t? | 
b Reference T5| 
c Reference "M \ 
^Reference 32 1 
e Reference 22 1 
-^Reference 3T| 
g Reference[30| 



TABLE II: Equilibrium distance, in A, of Ar on graphite and Al(100) and H2 on Al(100), computed 
using the standard DFT-revPBE calculation, and including the VdW correction, compared to the 
LDA result, and available theoretical and experimental (in parenthesis) reference data. 



system 


revPBE 


revPBE+VdW 


LDA 


ref. 


Ar-gr. hollow 


4.69 


3.61 


3.16 


3.33 a , 3.32 6 


Ar-gr. top 


4.95 


3.61 


3.22 


3.37 6 


Ar-gr. bridge 


4.69 


3.55 


3.22 


3.37 6 


Ar-benzene 


4.79 


3.57 


3.27 


3.41 c (3.68 d ) 


Ar-Al(lOO) hollow 


5.34 


4.66 


3.48 




Ar-Al(lOO) top 


5.34 


4.88 


3.57 




H 2 -A1(100) 


5.08 


4.62 


3.23 





"Reference 17]. 
b Reference 16]. 
c Referencc 32] . 
d Reference[33|. 




distance (A) 



